This script creates a merged list of treatments from both CalMapper and PFIRS with: - duplicates removed from PFIRS, defined as within 1 km and 30 days from CalMapper record - years 2020-2022 - filtered to AEPP threshold - buffered overlay method to account for PFIRS being point data
file = "Rdata/CalMapper_activities_fire.Rdata"
load(file)
file = "Rdata/PFIRS_2017-2022_pull2023.Rdata"
load(file)
cm %>%
group_by(YEAR) %>%
ggplot(aes(x = as.factor(YEAR), y = ACTIVITY_START))+
geom_point(position = position_jitter())
cm %>%
group_by(YEAR) %>%
ggplot(aes(x = as.factor(YEAR), y = ACTIVITY_END))+
geom_point(position = position_jitter())
Yup, it doesn’t look like there are issues
cm <- cm %>%
filter(YEAR %in% seq(2020, 2022))
pfirs %>%
group_by(Year) %>%
ggplot(aes(x = as.factor(Year), y = Burn_Date))+
geom_point(position = position_jitter())
pfirs <- pfirs %>%
filter(Year %in% seq(2020, 2022))
cm %>%
st_drop_geometry() %>%
group_by(ACTIVITY_DESCRIPTION) %>%
count()
## # A tibble: 3 × 2
## # Groups: ACTIVITY_DESCRIPTION [3]
## ACTIVITY_DESCRIPTION n
## <chr> <int>
## 1 Broadcast Burn 626
## 2 Cultural Burning 4
## 3 Pile Burning 1445
nrow_old <- nrow(cm)
cm <- cm %>%
filter((ACTIVITY_DESCRIPTION == "Pile Burning" & TREATED_ACRES >= 25) | (ACTIVITY_DESCRIPTION %in% c("Broadcast Burn", "Cultural Burning") & TREATED_ACRES >= 50))
deleted_n <- nrow_old - nrow(cm)
print(paste("Deleted", deleted_n, "records below the size threshold"))
## [1] "Deleted 1689 records below the size threshold"
pfirs %>%
st_drop_geometry() %>%
group_by(Burn_Type) %>%
count()
## # A tibble: 11 × 2
## # Groups: Burn_Type [11]
## Burn_Type n
## <chr> <int>
## 1 Broadcast 1526
## 2 Hand Pile 3087
## 3 Hand Piles 43
## 4 Landing Pile 226
## 5 Landing Piles 27
## 6 Machine Pile 1484
## 7 Machine Piles 2
## 8 Multiple Fuel Types 61
## 9 Multiple Fuels 46
## 10 UNK 41
## 11 Unknown 329
pfirs <- pfirs %>%
mutate(Burn_Type_Simple = case_when(
Burn_Type %in% c("Hand Pile", "Hand Piles", "Landing Pile", "Landing Piles", "Machine Pile", "Machine Piles") ~ "Pile",
Burn_Type == "Broadcast" ~ "Broadcast",
Burn_Type %in% c("Multiple Fuel Types", "Multiple Fuels", "UNK", "Unknown") ~ "Unknown/Multiple",
TRUE ~ NA
))
pfirs %>%
st_drop_geometry() %>%
group_by(Burn_Type, Burn_Type_Simple) %>%
count()
## # A tibble: 11 × 3
## # Groups: Burn_Type, Burn_Type_Simple [11]
## Burn_Type Burn_Type_Simple n
## <chr> <chr> <int>
## 1 Broadcast Broadcast 1526
## 2 Hand Pile Pile 3087
## 3 Hand Piles Pile 43
## 4 Landing Pile Pile 226
## 5 Landing Piles Pile 27
## 6 Machine Pile Pile 1484
## 7 Machine Piles Pile 2
## 8 Multiple Fuel Types Unknown/Multiple 61
## 9 Multiple Fuels Unknown/Multiple 46
## 10 UNK Unknown/Multiple 41
## 11 Unknown Unknown/Multiple 329
nrow_old <- nrow(pfirs)
pfirs <- pfirs %>%
filter((Burn_Type_Simple == "Pile" & Acres_Burned >= 25) | (Burn_Type_Simple %in% c("Broadcast", "Unknown/Multiple") & Acres_Burned >= 50))
deleted_n <- nrow_old - nrow(pfirs)
print(paste("Deleted", deleted_n, "records below the size threshold"))
## [1] "Deleted 5617 records below the size threshold"
if(object.size(cm) > object.size(pfirs)){
print("CalMapper is more GB than PFIRS")
}
## [1] "CalMapper is more GB than PFIRS"
pfirs <- st_transform(pfirs, st_crs(cm))
pfirs <- pfirs %>%
mutate(ID_CT = seq(1, nrow(pfirs)))
pfirs <- pfirs %>%
mutate(Data_source = "PFIRS")
cm <- cm %>%
mutate(Data_source = "CalMapper")
pfirs$Year <- as.factor(pfirs$Year)
cm$Year <- as.factor(cm$YEAR)
cm <- cm %>%
select(-YEAR)
pfirs_buff <- st_buffer(pfirs, dist = 1000)
st_write(pfirs_buff, dsn = "shapefiles", layer = "pfirs_buffer", driver = "ESRI Shapefile", delete_layer = T)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `pfirs_buffer' using driver `ESRI Shapefile'
## Writing layer `pfirs_buffer' to data source `shapefiles' using driver `ESRI Shapefile'
## Updating existing layer pfirs_buffer
## Writing 1255 features with 13 fields and geometry type Polygon.
st_crs(pfirs) == st_crs(cm)
## [1] TRUE
mapview(pfirs, zcol = "Year", col.regions = colors, alpha.regions = 0.5)+
mapview(cm, zcol = "Year", col.regions = colors, alpha.regions = 0.5)+
mapview(pfirs_buff, zcol = "Year", col.regions = colors, alpha.regions = 0.3)